home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cagdswep.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  12.1 KB  |  352 lines
  1. /******************************************************************************
  2. * CagdSwep.c - Sweep srf operator out of given cross section and an axis.     *
  3. *******************************************************************************
  4. * Written by Gershon Elber, May. 91.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. static void TransformCrossSection(CagdRType **SPoints, int Index,
  10.               CagdCrvStruct *CrossSection, CagdRType *Position,
  11.               CagdRType Scale, CagdVecStruct *Tangent,
  12.               CagdVecStruct *Normal);
  13. static void GenTransformMatrix(CagdMType Mat, CagdRType *Trans,
  14.                    CagdVecStruct *Normal, CagdVecStruct *Dir,
  15.                    CagdRType Scale);
  16.  
  17. /******************************************************************************
  18. * Constructs a sweep surface using the following curves:                      *
  19. * 1. CrossSection - defines the basic cross section of the sweep. Must be in  *
  20. *    the XY plane.                                  *
  21. * 2. Axis - a 3D curve the CrossSection will be swept along such that the     *
  22. *    Axis normal aligns with the Y axis of the cross section. If Axis is      *
  23. *    linear (i.e. no normal), the normal is picked randomly or to fit the non *
  24. *    linear part of the Axis (if any).                          *
  25. * 3. Scale - a scaling curve for the sweep, If NULL a scale of Scale is used. *
  26. * 4. Frame - a curve or a vector that specifies the orientation of the sweep  *
  27. *    by specifying the axes curve's normal. If Frame is a vector, it is a     *
  28. *    constant normal. If Frame is a curve (FrameIsCrv = TRUE), it is assumed  *
  29. *    to be a vector field normal. If NULL, it is computed from axes curve's   *
  30. *    pseudo Frenet frame, that minimizes rotation.                  *
  31. ******************************************************************************/
  32. CagdSrfStruct *CagdSweepSrf(CagdCrvStruct *CrossSection, CagdCrvStruct *Axis,
  33.                 CagdCrvStruct *ScalingCrv, CagdRType Scale,
  34.                 VoidPtr Frame, CagdBType FrameIsCrv)
  35. {
  36.     CagdSrfStruct
  37.     *Srf = NULL;
  38.     CagdVecStruct Normal, *Vec, TVec;
  39.     CagdPointType
  40.     SrfPType = CAGD_PT_E3_TYPE;
  41.     CagdGeomType
  42.     SrfGType = CAGD_SBSPLINE_TYPE;
  43.     int i, j, k,
  44.     ULength = CrossSection -> Length,
  45.     VLength = Axis -> Length,
  46.     UOrder = CrossSection -> Order,
  47.     VOrder = Axis -> Order;
  48.     CagdRType **Points,
  49.     *AxisKV = NULL,
  50.     *AxisNodes = NULL,
  51.     *AxisNodePtr = NULL,
  52.     *AxisWeights = NULL,
  53.     *FrameVec = FrameIsCrv ? NULL : (CagdRType *) Frame;
  54.     CagdCrvStruct
  55.     *FrameCrv = FrameIsCrv ? (CagdCrvStruct *) Frame : NULL;
  56.  
  57.     switch (Axis -> GType) {
  58.     case CAGD_CBEZIER_TYPE:
  59.         SrfGType = CAGD_SBEZIER_TYPE;
  60.         AxisKV = BspKnotUniformOpen(VLength, VOrder, NULL);
  61.         AxisNodes = AxisNodePtr = BspKnotNodes(AxisKV,
  62.                            VLength + VOrder,
  63.                            VOrder);
  64.         break;
  65.     case CAGD_CBSPLINE_TYPE:
  66.         SrfGType = CAGD_SBSPLINE_TYPE;
  67.         AxisKV = Axis -> KnotVector;
  68.         AxisNodePtr = AxisNodes = BspKnotNodes(Axis -> KnotVector,
  69.                            VLength + VOrder,
  70.                            VOrder);
  71.         break;
  72.     case CAGD_CPOWER_TYPE:
  73.         FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
  74.         break;
  75.     default:
  76.         FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
  77.         break;
  78.     }
  79.     if (CAGD_IS_RATIONAL_CRV(Axis))
  80.     AxisWeights = Axis -> Points[W];
  81.  
  82.     switch (CrossSection -> GType) {
  83.     case CAGD_CBEZIER_TYPE:
  84.         break;
  85.     case CAGD_CBSPLINE_TYPE:
  86.         SrfGType = CAGD_SBSPLINE_TYPE;
  87.         break;
  88.     case CAGD_CPOWER_TYPE:
  89.         FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
  90.         break;
  91.     default:
  92.         FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
  93.         break;
  94.     }
  95.  
  96.     if (ScalingCrv) {
  97.     int ScaleKVLen,
  98.         AxisKVLen = Axis -> Order + Axis -> Length;
  99.  
  100.     switch (ScalingCrv -> GType) {
  101.         case CAGD_CBEZIER_TYPE:
  102.         ScalingCrv = CnvrtBezier2BsplineCrv(ScalingCrv);
  103.             break;
  104.         case CAGD_CBSPLINE_TYPE:
  105.             ScalingCrv = CagdCrvCopy(ScalingCrv);
  106.         break;
  107.         case CAGD_CPOWER_TYPE:
  108.         FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
  109.         break;
  110.         default:
  111.         FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
  112.         break;
  113.     }
  114.  
  115.     ScaleKVLen = ScalingCrv -> Order + ScalingCrv -> Length;
  116.  
  117.         /* Affine trans. ScalingCrv KV to match Axis. */
  118.     BspKnotAffineTrans(ScalingCrv -> KnotVector, ScaleKVLen,
  119.                AxisKV[0] - ScalingCrv -> KnotVector[0],
  120.                (AxisKV[AxisKVLen - 1] - AxisKV[0]) /
  121.                (ScalingCrv -> KnotVector[ScaleKVLen - 1] -
  122.                 ScalingCrv -> KnotVector[0]));
  123.     }
  124.  
  125.     if (FrameCrv) {
  126.     int FrameKVLen,
  127.         AxisKVLen = Axis -> Order + Axis -> Length;
  128.  
  129.     switch (FrameCrv -> GType) {
  130.         case CAGD_CBEZIER_TYPE:
  131.         FrameCrv = CnvrtBezier2BsplineCrv(FrameCrv);
  132.             break;
  133.         case CAGD_CBSPLINE_TYPE:
  134.             FrameCrv = CagdCrvCopy(FrameCrv);
  135.         break;
  136.         case CAGD_CPOWER_TYPE:
  137.         FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
  138.         break;
  139.         default:
  140.         FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
  141.         break;
  142.     }
  143.  
  144.     FrameKVLen = FrameCrv -> Order + FrameCrv -> Length;
  145.  
  146.         /* Affine trans. FrameCrv KV to match Axis. */
  147.     BspKnotAffineTrans(FrameCrv -> KnotVector, FrameKVLen,
  148.                AxisKV[0] - FrameCrv -> KnotVector[0],
  149.                (AxisKV[AxisKVLen - 1] - AxisKV[0]) /
  150.                (FrameCrv -> KnotVector[FrameKVLen - 1] -
  151.                 FrameCrv -> KnotVector[0]));
  152.     }
  153.  
  154.     if (CAGD_IS_RATIONAL_CRV(Axis) || CAGD_IS_RATIONAL_CRV(CrossSection))
  155.     SrfPType = CAGD_PT_P3_TYPE;
  156.  
  157.     switch (SrfGType) {
  158.     case CAGD_SBEZIER_TYPE:
  159.         Srf = BzrSrfNew(ULength, VLength, SrfPType);
  160.         break;
  161.     case CAGD_SBSPLINE_TYPE:
  162.         Srf = BspSrfNew(ULength, VLength, UOrder, VOrder, SrfPType);
  163.         if (CrossSection -> GType == CAGD_CBSPLINE_TYPE)
  164.         CAGD_GEN_COPY(Srf -> UKnotVector, CrossSection -> KnotVector,
  165.                   sizeof(CagdRType) * (ULength + UOrder));
  166.         else
  167.         BspKnotUniformOpen(ULength, UOrder, Srf -> UKnotVector);
  168.         if (Axis -> GType == CAGD_CBSPLINE_TYPE)
  169.         CAGD_GEN_COPY(Srf -> VKnotVector, Axis -> KnotVector,
  170.                   sizeof(CagdRType) * (VLength + VOrder));
  171.         else
  172.         BspKnotUniformOpen(VLength, VOrder, Srf -> VKnotVector);
  173.         break;
  174.     default:
  175.         FATAL_ERROR(CAGD_ERR_UNDEF_SRF);
  176.         break;
  177.     }
  178.     Points = Srf -> Points;
  179.  
  180.     /* Compute the first normal to be used to orient first cross section. */
  181.     if (FrameVec == NULL && FrameCrv == NULL) {
  182.     /* Compute the normal to axis curve and use it to align the +Y axis  */
  183.     /* of cross section with that vector. If Axis curve has no normal    */
  184.     /* (i.e. it is a linear segment), an arbitrary normal is picked.     */
  185.     Vec = CagdCrvNormal(Axis, *AxisNodePtr);
  186.     if (Vec != NULL) {
  187.         CAGD_COPY_VECTOR(Normal, *Vec);
  188.     }
  189.     else {
  190.         Vec = CagdCrvTangent(Axis, *AxisNodePtr);
  191.         Normal.Vec[0] = Normal.Vec[1] = Normal.Vec[2] = 0.0;
  192.         if (ABS(Vec -> Vec[0]) <= ABS(Vec -> Vec[1]) &&
  193.         ABS(Vec -> Vec[0]) <= ABS(Vec -> Vec[2]))
  194.         Normal.Vec[0] = 1.0;
  195.         else if (ABS(Vec -> Vec[1]) <= ABS(Vec -> Vec[0]) &&
  196.              ABS(Vec -> Vec[1]) <= ABS(Vec -> Vec[2]))
  197.         Normal.Vec[1] = 1.0;
  198.         else
  199.         Normal.Vec[2] = 1.0;
  200.  
  201.         CROSS_PROD(TVec.Vec, Vec -> Vec, Normal.Vec);
  202.         CAGD_NORMALIZE_VECTOR(TVec);
  203.         CROSS_PROD(Normal.Vec, TVec.Vec, Vec -> Vec);
  204.     }
  205.     }
  206.  
  207.     /* For each ctl points of the axis, transform the cross section          */
  208.     /* according to ctl point position, tangent to axis at the point and in  */
  209.     /* such a way to minimize Normal change.                     */
  210.     for (i = 0; i < VLength; i++, AxisNodePtr++) {
  211.     CagdRType PosE3[3], ScaleE2[2],
  212.         *Scaling = ScalingCrv ? CagdCrvEval(ScalingCrv, *AxisNodePtr) : NULL;
  213.     CagdVecStruct
  214.         *Tangent = CagdCrvTangent(Axis, *AxisNodePtr);
  215.  
  216.     if (Scaling)
  217.         CagdCoerceToE2(ScaleE2, &Scaling, -1, ScalingCrv -> PType);
  218.     else
  219.         ScaleE2[1] = Scale;
  220.     if (AxisWeights)
  221.         ScaleE2[1] /= AxisWeights[i];
  222.  
  223.     /* If Normal is fully specified, get it now. */
  224.     if (FrameVec != NULL) {
  225.         PT_COPY(Normal.Vec, FrameVec);
  226.     }
  227.     else if (FrameCrv != NULL) {
  228.         CagdRType
  229.         *FrameCrvVal = CagdCrvEval(FrameCrv, *AxisNodePtr);
  230.  
  231.         CagdCoerceToE3(Normal.Vec, &FrameCrvVal, -1, FrameCrv -> PType);
  232.         CAGD_NORMALIZE_VECTOR(Normal);
  233.     }
  234.  
  235.     CagdCoerceToE3(PosE3, Axis -> Points, i, Axis -> PType);
  236.     TransformCrossSection(Points, i * ULength, CrossSection,
  237.                   PosE3, ScaleE2[1], Tangent, &Normal);
  238.     }
  239.  
  240.     /* Do fixups if axis is a rational curve (note surface is P3). */
  241.     if (AxisWeights) {
  242.     if (CAGD_IS_RATIONAL_CRV(CrossSection)) {
  243.         /* Need only scale by the Axis curve weights: */
  244.         for (j = 0, k = 0; j < VLength; j++)
  245.         for (i = 0; i < ULength; i++, k++) {
  246.             Points[X][k] *= AxisWeights[j];
  247.             Points[Y][k] *= AxisWeights[j];
  248.             Points[Z][k] *= AxisWeights[j];
  249.             Points[W][k] *= AxisWeights[j];
  250.         }
  251.     }
  252.     else {
  253.         /* Weights do not exists at the moment - need to copy them. */
  254.         for (j = 0, k = 0; j < VLength; j++)
  255.         for (i = 0; i < ULength; i++, k++) {
  256.             Points[X][k] *= AxisWeights[j];
  257.             Points[Y][k] *= AxisWeights[j];
  258.             Points[Z][k] *= AxisWeights[j];
  259.             Points[W][k] = AxisWeights[j];
  260.         }
  261.     }
  262.     }
  263.  
  264.     if (Axis -> GType == CAGD_CBEZIER_TYPE)
  265.         IritFree((VoidPtr) AxisKV);
  266.     if (ScalingCrv)
  267.     CagdCrvFree(ScalingCrv);
  268.     IritFree((VoidPtr) AxisNodes);
  269.  
  270.     return Srf;
  271. }
  272.  
  273. /******************************************************************************
  274. * Transform The CrossSection Points, to Position such that Tangent is         *
  275. * perpendicular to the cross section (which is assumed to be on the XY        *
  276. * plane). The +Y axis of the cross section is aligned with Normal direction   *
  277. * to minimize twist along the sweep and been updated to new normal.          *
  278. *   Transformed cross section is place into srf Points, SPoints starting from *
  279. * index SIndex.                                      *
  280. *   All agrument vectors are assumed to be normalized to a unit length.          *
  281. ******************************************************************************/
  282. static void TransformCrossSection(CagdRType **SPoints, int SIndex,
  283.               CagdCrvStruct *CrossSection, CagdRType *Position,
  284.               CagdRType Scale, CagdVecStruct *Tangent,
  285.               CagdVecStruct *Normal)
  286. {
  287.     CagdPointType
  288.     PType = CrossSection -> PType;
  289.     CagdBType
  290.     IsNotRational = !CAGD_IS_RATIONAL_PT(PType);
  291.     CagdVecStruct TVec;
  292.     CagdMType Mat;
  293.     CagdCrvStruct
  294.     *CSCopy = CagdCrvCopy(CrossSection);
  295.     int i, j,
  296.     MaxCoord = CAGD_NUM_OF_PT_COORD(PType),
  297.     Len = CSCopy -> Length;
  298.     CagdRType R,
  299.     **CSPoints = CSCopy -> Points;
  300.  
  301.     /* Fix the Normal to be perpendicular to the Tangent vector is a minimal */
  302.     /* rotation. This ensures minimal twist in the resulting surface.         */
  303.     R = DOT_PROD(Normal -> Vec, Tangent -> Vec);
  304.     TVec = *Tangent;
  305.     CAGD_MULT_VECTOR(TVec, R);
  306.     CAGD_SUB_VECTOR(*Normal, TVec);
  307.     CAGD_NORMALIZE_VECTOR(*Normal);
  308.  
  309.     GenTransformMatrix(Mat, Position, Normal, Tangent, Scale);
  310.     CagdCrvMatTransform(CSCopy, Mat);
  311.  
  312.     for (i = 0; i < Len; i++)
  313.     for (j = IsNotRational; j <= MaxCoord; j++)
  314.         SPoints[j][SIndex + i] = CSPoints[j][i];
  315.  
  316.     CagdCrvFree(CSCopy);
  317. }
  318.  
  319. /******************************************************************************
  320. *   Routine to preper transformation martix to do the following (in this      *
  321. * order): scale by Scale, rotate such that the Z axis is in direction Dir     *
  322. * and Y is colinear with the Normal and then translate by Trans.          *
  323. *    Algorithm: given the Trans vector, it forms the 4th line of Mat. Dir is  *
  324. * used to form the second line (the first 3 lines set the rotation), and      *
  325. * finally Scale is used to scale first 3 lines/columns to the needed scale:   *
  326. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord  *
  327. *                |  Bx  By  Bz  0 |  system into T, N & B as required and     *
  328. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |  then translate it to C. T, N, B are      *
  329. *                |  Cx  Cy  Cz  1 |  scaled by Scale.                  *
  330. * N is exactly Dir (unit vec). T is set to be the Normal and B their cross    *
  331. * product.                                      *
  332. *   All agrument vectors are assumed to be normalized to a unit length.          *
  333. ******************************************************************************/
  334. static void GenTransformMatrix(CagdMType Mat, CagdRType *Trans,
  335.                    CagdVecStruct *Normal, CagdVecStruct *Dir,
  336.                    CagdRType Scale)
  337. {
  338.     int i;
  339.     CagdVecStruct B;
  340.  
  341.     CROSS_PROD(B.Vec, Dir -> Vec, Normal -> Vec);
  342.  
  343.     for (i = 0; i < 3; i++) {
  344.     Mat[0][i] = Normal -> Vec[i] * Scale;
  345.     Mat[1][i] = B.Vec[i] * Scale;
  346.     Mat[2][i] = Dir -> Vec[i] * Scale;
  347.     Mat[3][i] = Trans[i];
  348.     Mat[i][3] = 0.0;
  349.     }
  350.     Mat[3][3] = 1.0;
  351. }
  352.